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CONSTRUCTING ELLIPTIC CURVES OVER FINITE FIELDS 
WITH PRESCRIBED TORSION 

ANDREW V. SUTHERLAND 



Abstract. The modular curve Xi{N) parametrizes elliptic curves with a 
point of order N. For Af < 50 we obtain plane models of Xi{N) that have 
been optimized for fast computation, and provide explicit birational maps to 
transform a point on our model of Xi (N) to an elliptic curve. Over a finite 
field, these allow us to quickly construct elliptic curves containing a point of 
order N, and can accelerate the search for an elliptic curve whose order is 
divisible by N. 



1. Introduction 

By Mazur's theorem [12 , the order of a nontrivial torsion point on an eUiptic 
curve over the rational numbers must belong to the set 

r = {2,3,4,5,6,7,8,9,10, 12}. 

Conversely, for each e T, an infinite family of elliptic curves over Q containing 
a point of order N is exhibited by the parametrizations of Kubert P3. Q Over a 
finite field F^, these parametrizations provide an easy way to generate universal 
families of curves whose order is divisible by N. This can accelerate applications 
that search for an elliptic curve with a particular property, such as a curves with 
smooth order (as in the elliptic curve factorization method [1]), or curves with a 
particular endomorphism ring (as when computing Hilbert class polynomials with 
the Chinese Remainder Theorem 2 ). 

To generate an elliptic curve E/¥q with non-trivial 7-torsion, for example, one 
applies [TT]. Pick r € F^, then use h — — and c = r'^ — r to define 

(1) E{b, c) : + (1 - c)xy -by^x^~ bx'^ . 

Provided E{b, c) is nonsingular, we obtain an elliptic curve on which the point 
P = (0, 0) has order 7. By contrast, obtaining such a curve by trial and error is far 
more time consuming: testing for 7-torsion typically involves finding the roots of a 
degree 24 polynomial (the 7-division polynomial), and several curves may need to 
be tested (approximately six on average) . 

Mazur's theorem limits us to iV € T, but we can proceed further if we do not re- 
strict ourselves to curves defined over Q. Reichert treats N G {11, 13, 14, 15, 16, 18} 
over quadratic extensions of Q using Xi (N) , the modular curve that parametrizes 
elliptic curves with a point of order TV [TS]. We may be able to realize a curve 
defined over K = Q[Vd] in F^, but only if d is a quadratic residue in Fq. 



2000 Mathematics Subject Classification. Primary 14H52; Secondary 11G20. 
^Kubert also addresses the torsion subgroups Z/2Z X Z/2AZ for N = 1,2,3,4. We consider 
only the subgroups Z/NZ here. 
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Alternatively, we can use a point on Xi{N)/¥q to directly construct E{b,c)/Wq 
containing a point of order N. This applies in any sufficiently large finite field, for 
all N > 3 (see [TT] for < 3). For A e T the curve Xi{N) has genus and 
we simply obtain the Kubert parametrizations, but in general we construct E{b, c) 
from a point {x,y) on Ai(A^) via a birational map that depends on the defining 
equation we choose for Ai(A^). 

For example, to obtain a curve with non-trivial 13-torsion, we use a point on 

which may be obtained by choosing x Cz¥q at random and attempting to solve the 
resulting quadratic equation for y in We then apply the transformation 

r — I — xy, 

s = l~ xy/{y + l), 

set c = s(r — 1) and b = cr, and construct E{b, c). If we obtain a singular curve (or 
if y = — 1) we try again with a different point on Ai(I3) (this rarely happens). 

To apply this method we require a defining equation for Xi{N)/¥q along with 
a suitable birational map. For fast computation we seek a plane model f{x, y) = 
that minimizes the degree d of one of its variables. For A^ < 18 one can derive these 
from the results of Reichert (and Kubert). Reichert's method can be applied to 
N > 18, but the "raw" form of Ai(A^) initially obtained is quite large and of higher 
degree than necessary. More compact defining equations for Ai(A^) are given by 
Yang |21j for A^ < 22, but these do not minimize d. The minimal value d — d{N) 
is a topic of some interest [U [7l |9j [131 [13 1 since we can construct (infinitely many) 
elliptic curves containing a point of order N over number fields of degree d. For 
A^ > 18, few values of d{N) are known (see sequence A146879 in the OEIS [T7]). 

Given a plane model for Ai(A), we may attempt to reduce its complexity (de- 
gree, number of terms, and coefficient size) through a judiciously chosen sequence 
of rational transformations. This procedure is somewhat ad hoc, however, and find- 
ing an optimal (or even good) sequence becomes difficult for larger A^. We treat 
this as a combinatorial optimization problem, applying standard search techniques 
to obtain a solution that is locally optimal under a relation we define. We cannot 
claim that the results are globally optimal, but they do yield an upper bound on 
d{n) . For N < 22 we are able to match known lower bounds for d{N) O |7l [9] , 
including d{l9) = 5, which we believe to be newH Results for A^ < 30 are listed in 
the appendix, and are available in electronic form for A^ < 50. 

For odd A^ we also show how to efficiently generate E/¥q with a point of order 
AN, or satisfying ^E = 2N mod AN, using our results for Ai(2A^). 

2. Computing the raw form of Ai(A^) 

Following Reichert [15] . we summarize the method to obtain a plane model for 
Xi{N) in the form F{r,s) = 00 The equation Eib,c) in ^ is the Tate normal 
form of an elliptic curve (called a Kubert curve in [1]). Any elliptic curve containing 



When Xi{N) has genus 1 (A'^ = 11, 14, 15) we may obtain additional points more efiiciently 
using the group operation on Xi{N) (see Section |4)l. 
■^But wc do not achieve d{24) = 4 imphed by 6 . 

^Reichert uses auxiUary variables m = s(l — r)/(l — s) and t = (r — s)/(l — s). We find it 
preferable to work directly with r and s. 



CONSTRUCTING ELLIPTIC CURVES WITH PRESCRIBED TORSION 



3 



a point of order greater than 3 can be put in this form (see V.5 of [10]). The 
discriminant of E{b, c) is 

(2) A{b, c) = b^{16b^ - 8bc^ - 206c + b + c{c- if). 

To ensure E(b, c) is nonsingular we require A(6, c) 7^ 0, so assume b ^ 0. Applying 
the group law for elliptic curves [TBI 111.2.3], we double the point P — (0,0) to 
obtain 2P = {b,bc), and for n > 1 compute the point (n + 1)P — (.Tn+i, J/n+i) hi 
terms of nP = (a;„, y„) using 

(3) Xn+l = byn/xl, Vn+l = 6^ (x^ - yn)/xl. 

We find that the inverse of nP — (a;„, ?;„) is 

(4) -nP = {x„,b+ic-l)xn~yn)■ 
l^ P is an A^-torsion point and m + ji = N, then we must have mP = —nP. If 
m ^ n this implies Xm — Xn, and if to = n we have 2j/„ = 6 + (c — l)a;„. For 6^0 
this requires N > 3. When Xm = x„ either mP = nP or mP = —nP, and in the 
latter case P is an A-torsion point. If we choose to = f-^^^^] and n — [-^^^^J we 
ensure that mP ^ nP, obtaining a necessary and sufficient condition for A-torsion: 

(5) NP = Oe ^ x„, = Xn, 

valid for A^ > 3. The first three multiples of P are: 

P-(0,0), 
2F=(fe,6c), 
3F= (c,5-c). 

We see immediately that P is a point of order 4 exactly when c = 0, and P is a 
point of order 5 exactly when b = c. For N > 5 define: 

r — b/c, b — rs{r — 1), 

s ~ (? I {b — c), c = s(r — 1), 

and note that r ^ {0, 1}, and s ^ 0. 

We now apply ([3]) to compute x„ in terms of r and s. Values for n < 10 are 
listed in Table [TJ To obtain the raw form of Ai ( A^) , we start with the equation 
Xm = Xn from ([5]), then clear denominators and subtract to obtain an equation of 
the form F*{r,s) = 0, where the polynomial F*{r,s) has integer coefficients. We 
then remove from F* any factors prohibited by our assumptions (namely r, r — 1, 
and s), and also factors corresponding to M-torsion for any M > 5 dividing A^ 
(such as s — 1 for A/ = 6 and r — s for A/ = 7) 

Let F{r, s) denote the polynomial that remains. We claim that F{r, s) = is 
a defining equation for Ai(A^). By construction, any solution to F{r,s) — will 
produce a curve E{b, c), with c — s{r— 1) and b = rc, on which P is a point of order 
A^, provided that A(fe, c) 7^ 0. Conversely, any curve E{b,c) on which P has order 
A^ > 5 yields a solution r = b/c, s = c^/{b — c) to F{r,s) = 0. These statements 
hold for any field K, provided that we verify A (6, c) 7^ in K. 



More generally, these can be recognized by computing the raw form of each Xi(M). In 
practice F(r, s) is simply the largest irreducible factor of F*{r, s). 
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X2 = rs{r — 1) 
X3, = s{r - 1) 
Xi = r{r — 1) 
0:5 = rs(s — 1) 

xq = s{r — l){r ~ s) I (s — 1)^ 

X7 = rs(r — l)(s — l)(rs — 2r + 1) / (r — s)^ 

xg = r(r - l)(r - s)(r - + s - 1) / (rs - 2r + 1)^ 

a;g = s{r — l)(rs — 2r + l)(rs^ — 3rs + r + s^) / (r — + s — 1)^ 

a^io = — + s — l)(r^ — rs^ + Srs^ — 4rs + s) / (rs^ — 3rs + r + s^)^ 

Table 1. x-coordinates of nP for n < 10. 

When N — 16, for example, putting xg = xi in the form F*{r, s) — yields 

F*{r,s) ^s{r - l)(r - s)^(rs - 2r + l)(rs^ - 3rs + r + s + s^) 

- rs{r - l)(s - l)(rs - 2r + l)(r - + s - 1)^. 

The nonzero factors s and r — 1 may be removed, and also the factor rs — 2r + 1, 
which can be zero only when P has order 8. Thus we obtain 

F{r, s) = {r ~ sY{rs^ - 3rs + r + s + s^) - r{s - l)(r - + s - 1)^. 

When expanded, this yields the entry for N — 16 in Table |4l The polynomials 
F{r, s) for N up to 50 are available in electronic form from the author. The largest 
of these has 1,791 terms and maximum coefficient on the order of 10^^. 

3. Reducing the complexity of F{r, s) = 

To facilitate fast computation we wish to simplify the raw from of Xi{N). We 
seek a birationally equivalent curve f{x, y) = that minimizes the degree of one of 
its variables (say y). Subject to this constraint, we would like to make / monic in 
y and also to minimize the degree in x, the number of terms, and the size of the 
coefficients (roughly in that order of priority). One typically approaches this prob- 
lem by attempting to remove singularities from i^(r, s) through a combination of 
translations and inversions (see |15j for examples). We take a more naive approach 
that allows us to easily automate the process. 

There are three basic types of transformations we will use: 

(1) Translate: x x + a or y y + a. 

(2) Invert: x 1/x or y 1/y. 

(3) Separate: x--~*l/x, y-^yjx or x-^xjy, y~^l/y. 

These are clearly all invertible operations. The third type combines an inversion 
and a division, but we find it works well as an atomic unit. In order to bound the 
number of atomic operations, we let a € {±1}, giving a total of eight. 

Consider the directed graph G on the set C of plane curves that can be obtained 
from F(r, s) = by applying a finite sequence of the transformations above, with 
edges labeled by the corresponding operation. A path in G defines a birational map 
(the composition of the operations labeling its edges), and any path can be reversed 
to yield the inverse map. Starting from the curve Co defined by F(r, s) = 0, we 
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want to find a path to a "better" curve C . To make this precise, we associate to 
each integer polynomial f{x,y) a vector of nonnegative integers 

«(/) = {dy,my,dx,dtot,t,S), 

whose components are defined by: 

• dy is the degree of / in y; 

• niy is if no term of / is a multiple of xy'^^ and 1 otherwise; 

• dtot is the total degree of /; 

• t is the number of terms in /; 

• S" is the sum of the absolute values of the coefficients of /. 

The component my will be zero exactly when / can be made monic as a polynomial 
in y. We order the vectors v{f) lexicographically, and to each C G C assign the vec- 
tor v{C) = miii{v{f{x,y)),v{f{y,x))}, where f{x,y) = defines C. We compare 
curves by comparing their vectors, obtaining a prewellordering of C. In particular, 
any subset of C contains a (not necessarily unique) minimal element. 

We now give a simple algorithm to search the graph G for a curve that is locally 
optimal within a radius R. We use -/V(C, k) to denote the set of curves connected 
to C by a path of length at most k in G. For G' € N{G, k) we let 0(C", G) denote 
the birational map defined by the path from G' back to G . 

1. Set C <— Co, fc = 1, and let be the identity map. 

2. While k<R: 

a. Determine a minimal element G' of N(C, fc). 

b. If v{G') < viG), then set ip ^ ip o (f>{G' , G), G ^ G', and fc ^ 0. 

c. Set fc ^ fc + 1. 

3. Output Ci = C and ip. 

The curve Ci output by the algorithm is our optimized plane model for Xi{N). It 
is birationally equivalent to the curve Gq defined by s) = 0, and the map ip 
carries points on Ci to points on Gq. 

To enumerate the neighbors of the curve G defined by f{x, y) = 0, the algorithm 
applies each of the eight atomic operations. The result of applying the birational 
map 4> with inverse ip is computed by expanding f{(j)xix,y),(/)y{x,y)) as a formal 
substitution of variables and clearing any denominators that result. Thus the trans- 
lation X a; — 1 is obtained by expanding f{x + 1, y), and the inversion x 1/x 
effectively replaces in f{x,y) with x''^"*. To enumerate N{G,k) involves ap- 
plying up to 8'^ possible sequences of operations (this number can be reduced by 
eliminating obviously redundant sequences), so the bound R cannot be very large. 
We have tested up to i? = 10, but find that R = 8 suffices to obtain the results 
given here. When R = 8 the algorithm takes less than an hour (on a 2.8 GHz AMD 
Athlon processor) for N < 50. 

Table [2] illustrates the algorithm's execution for N = 16. We begin with the 
curve Co defined by F{r, s) = 0, as listed in TablcHl and set C = Co with /(x, y) = 
F{x, y). The algorithm finds v{G) = v{f{y, x)) = (3, 1, 5, 6, 13, 40), indicating that 
the f{x,y) has degree 3 in x (in this case v{f{y,x)) is less than v{f{x,y)) so the 
roles of X and y are reversed). Additionally, /(x, y) is not monic in x, has degree 5 
in y, total degree 6, 13 terms, and the absolute values of its coefficients sum to 40. 

No curves within a distance fc = 1 are found that improve v{G), but for fc = 2 a 
curve C is found that is monic in x (and also degree 3), which implies v(C') < v{G). 
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Steps C : fix, y) ^ v{C) 

x^y'^ - Ax^y + 2x^ + 3x'^y'^ + 2x'^y -2x'^ - xy^ +Axy'^ (3,1,5,6,13,40) 

— lOxy^ + 6xy^ — 3xy + x + y'^ 

5,8 x^ + x^y^ - Sx'^y^ + 6x'^y^ - lOx^y^ + ix^y - x^ -2xy^ (3,0,7,7,13,40) 

+ 2xy^ + 3xy'^ + - 4y^ + j/^ 

1.3.8.6 x3 + a;2y4 + 2x2y3 + 4a;2y2 _ 52,2 _ 2a;t/'' - 8xy3 - I3a;y2 (3^0,4,6,13,68) 

+ 8a; + 2y^ + 8y^ + lOy^ - 4 

1 x^ +x^y'^ + 2x'^y^ + Ax'^y^ -2x^ - Axy-^ -5xy'^ + x (3,0,4,6,11,24) 

+ y^ + 2y3 + j/2 

1 x^ +x'^y* + 2x'^y^ + 4x'^y'^ + x'^ + 2xy-^ + 3xy'^ + 2y^ (3,0,4,6,8,16) 

5,6 2x^ + ix^y^ + 2x'^ + xy"" + Axy'^ + 2a;y + x + (3,0,4,5,8) 

2,4,5,6,8 -x^ + x^y^ -Ax'^y^ +Ax^y^2x'^ ^ixy^ ~%xy-x^2y (3,0,3,5,9,24) 

3 -a;3 + a;2y3 _ ^2y2 „ 2;2y + 3a;2 + 32:?/2 _ 4a; + 2y + 2 (3,3,0,5,9,18) 

4.5.1.7 -x^y^ -2x'^y~ x'^ ^xy^ + 2xy^ ^y"^ + Zy^ + 2y (2,1,3,4,8,13) 

4 -x^y^ + xy^ - xy^ ~ xy + x + y^ -y (2,1,3,4,7,7) 
8 -x^ + xy'^ -xy"^ -xy + x-y'^ + y (2,0,3,4,7,7) 
1 -x"^ + xy'^ - xy'^ - xy - x - y"^ (2,0,3,4,6,6) 



Table 2. Optimization of Xi (16). 

l:x^a; — 1, 2:a;^x + l, 3:?/^y— 1, A : y y + 1 
5 : x ~^ 1/x, 6:y--~*l/y, 7:x--~^l/x, y-^y/x, S:x-^x/y, y-^l/y. 

C is a minimal curve in N{C, 2), so C is replaced by C" and the map (p becomes 

x^y/x, y-^l/y. 

This reverses the sequence of steps 5,8 (as identified in the key to Table [2]) used 
to reach C from Co (so maps points on C" back to points on Co). The next 
improvement occurs when fc = 4. In this case reversing the path 1,3,8,6 from C to 
C yields the sequence 6,8,4,2, and ip becomes 

x-^ {y + l)l{xy + l), y^l/{y+\). 

The algorithm continues in this fashion, finding the sequence of curves listed in 
Table [2l until it is unable to find a better curve within the maximum search radius 
R. The resulting curve has minimal degree in x rather than y, so we swap variables 
(and adjust signs) to obtain the optimized curve 

(6) Xi(16) : y"^ + {x^ +x^ ~x + l)y + x^ = Q, 

which appears in Table El Corresponding changes to ip yield the birational map 

(7) r = l + (y+l)/(.Ty + y2), s ^ I + [y + I) / [xy ~ y"^), 

listed in Table [71 which carries points on the curve in ^ to points on Co. 

Table ini shows the improvement in the minimal degree diCi) and the number 
of terms t(Ci) obtained when the initial curve Co is transformed to the locally 
optimal curve Ci output by the algorithm. For comparison, we also list the genus 
of Xi(iV), obtained from sequence A029937 in the OEIS [H] (see Theorem 1.1 of 
[5] for a general formula). 
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The search procedure described above can be appHed to any plane curve defined 
over Q, but its effectiveness depends largely on finding singularities with small 
integer coordinates. Empirically, this works well with Xi{N), but other applica- 
tions may wish to modify the list of atomic operations to incorporate more general 
translations. Alternative search strategies, such as simulated annealing, may also 
be worth investigating. 

4. Application to finite fields 

We can use the optimized form of Xi{N) to efficiently generate elliptic curves 
containing a point of order N over the finite field F^, as described in the introduc- 
tion. Here we briefiy address a few topics relevant to practical implementation. We 
assume that Ci is defined by f{x,y) = 0, with dy <dx, and consider how we may 
use f{x,y) to efficiently generate a set of m elliptic curves over F^, each containing 
a point of order N . 

Except for a small set of points (those leading to singular curves and those for 
which ip is undefined), there is a one-to-one correspondence between points on Ci 
and nonsingular curves in Tate normal form on which the point P = (0, 0) has order 
N (see Section [5]). For large q, each possible j-invariant in Fg (and each twist) is 
represented by an approximately equal number of curves in Tate normal form. It 
follows that we can obtain a (nearly) uniform distribution over isomorphism classes 
of elliptic curves defined over F^ containing a point of order N , provided that we 
have a uniformly distributed sample of points on f{x,y) = 0. 

When d > 2 it is not a trivial task to efficiently generate a sample with uniform 
distribution. It is impractical to test random solutions to f{x,y) = 0, so instead 
we pick Xi e ¥q at random and compute the roots yij (if any) of the degree d 
polynomial hi{y) — f{xi,y) over F,. For each root j/^ of hi we include the point 
{xi,yij) in our set of m points. Assuming d this gives us an approximately 

uniform distribution (if we used only one root of hi this would not be true), but the 
points obtained are not all independent. In practice this does not pose a problem. 
At most d points share a common x value, and after mapping the points back to 
F{r, s) = and constructing -E(6, c) it is very difficult to discern any relationship 
among the curvesQ With this approach we expect to compute the roots of m 
polynomials hi{y), on average, in order to obtain m points on f{x,y) = 0. 

When Xi{N) has genus 1, the curve f{x, y) = is an elhptic curve, and we may 
use a more efficient approach: select one point at random, then compute multiples 
of it via the group operation. We can generate m random multiples using 0{\ogq + 
m log g/ log log g) group operations via standard multi-exponentiation techniques 
[22] . or we can compute multiples in an arithmetic sequence using just m + 0{\ogq) 
group operations. The latter approach does not generate independent points, but 
it is highly efficient: only 0(1) operations in Fp are required per point (assuming 
m ^ log q) . During this computation it is convenient to work with a model for 
Xi (N) in short Weierstrass form. These are provided in Table [3l along with the 
corresponding maps back to F{r, s) = 0. 

Having generated a set of m points on f{x,y) — 0, we apply the appropriate 
birational map to obtain points on i^(r, s) = 0. When doing so, we invert the 

^Alternatively, we could obtain a uniform independent distribution using as most one root of 
each hi, provided we discard it with a certain probability depending on the number of roots hi 
has. We do not regard this as a practical solution. 
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2/' = 

r = 
s = 


= 

(y^ 
1 + 


- 432a; + 8208 
- 108)/216 
(y- 108)/(6a; + 72) 
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= x^ 


- 675a; + 13662 






r = 


1 + 


(108a; - 36?; + 3564)/(3a;2 - xy 


- 342a; + 75y + 999) 




s = 


(6a; 


- 234)/(9a;-2/- 135) 
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y'-- 
t ^ 

r = 
s = 


= x^ 
{6x 
1 - 
1 - 


- 27a; + 8694 
-90)(18a; + 62;-918) 
t/(a;2y - 189a;2 + A2xy - 4050a; 
t/lx'^y - 81a;2 + 6a;y - 3402a; ~ 


~3y2 + 441y-1701) 
3?;2 + 981y- 35721) 



Table 3. Short Weierstrass form of Xi{N) with genus 1. 



denominators in paraUel, via the usual Montgomery trick [31 Alg. 11.15]. We then 
compute (6, c) pairs (using c = s{r — 1) and b = rc). In a field of characteristic not 
2 or 3, we may convert the curve E{b, c) to the short Weierstrass form: 

(8) y^ =x^ +Ax + B. 

Let d = c — 1 and e ^ <P — Ah. Through the admissible change of variables 

(9) a; = 36a;'-3e, y = 216y' + 108(dx' + 6), 
we find that 

A = 27(246d - e^), B = 54(e^ - 366de + 216&^), 

and (3d, —1086) is a point of order on = a;'^ + Ax + 

At some point during the process described above, we need to check that the 
discriminant A of each curve obtained is nonzero. This is most efficiently done at 
the end using A = — 4yl'^ — 27 B^. This may result in fewer than m curves being 
generated, but we can always obtain more points on Xi {N) if necessary. 

5. Prescribing 4-torsion 

For odd A^, we can use Xi{2N) to generate elliptic curves which contain a point 
of order 4A^ over in a manner that may be more efficient than using Xi(4A^). 
Alternatively, we can generate curves which contain a point of order 2A^ but do 
not contain a point of order 4A^. These results rely on efficiently computing the 
4-torsion of an elliptic curve using a known a point of order 2, which we obtain 
from the point P — (0, 0) of order 2A^. For odd A^, a curve with a point of order A^ 
has a point of order 4A^ if and only if it has a point of order 4. 

In fact, we only need the x-coordinate of NP, which can be computed as de- 
scribed in Section [2] (see Table [T] for A^ < 10). It will be convenient to work 
with the short Weierstrass form ([8]), so we assume that the point NP has been 
translated via ([9]) to the 2-torsion point (3 — (a;o,0) on the curve E defined by 
y2 ^ j(a,) ^ x^ + Ax + B. 



For N £ T, parametrizations which additionally provide a point with infinite order over Q 
are considered by Atkin and Morain in jlj. 
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Our strategy is to use the value xo to determine whether E contains a point 
of order 4 or not. In the best case this requires only a single test for quadratic 
residuacity in F^, and even in the worst case, a square root and two tests for 
quadratic residuacity sufhce. If the result is not as desired, we discard E and test 
another curve with a point of order 2N. On average we expect to test two curves. 
This is typically faster than either finding a point on Xi{4N), or using Xi{N) and 
computing 4-torsion without a known point of order 2. 

Lemma 1. If a ~ {u,v) and (3 — (a;o,0) are points on a nonsingular elliptic curve 
E defined by = f{^) = x^ + Ax + B over a field of characteristic not 2 then 

2a = f3 -^=^ {u- xof ^ f{xo), 

where f'{x) = ix"^ + A. 

Proof. Ii2a — (3 then the duplication formula for elliptic curves [1^ p. 59] implies 

_ - 2Au'^ - 8Bu + A^ 
" 4(^3 + Au + B) ■ 

Therefore u must satisfy 

u'* - 4x0^^ - 2Au'^ - {AAxo + 8B)u - ABxq + = 0. 
Since P = {xq, 0) G E, we have x^ + Axq + B = 0. Substituting for B yields 

- Axqu^ - 2Au'^ + {8x1 + 4:Axo)u + Ax^ + AAx^ + A^. 
We now set u — z + xq and rewrite this as 

{z^-i3xl + A)f^0. 

Therefore 

{u - xoY ^3x1 + A = f'ixo), 

as desired. Reversing the argument yields the converse, provided f{u) ^ 0. But 
if u is a root of /, then one can show that (u — xof' — f'{xo) implies D{f) — 0, 
contradicting the fact that E is nonsingular. □ 

There may be 1 or 3 points of order 2 on E. The x-coordinates of the other two 
(if they exist) are the roots xi and X2 of f{x)/{x — xo), which we can determine 
with the quadratic formula. We now give our main result for treating 4-torsion. 

Proposition 1. Let (xq, 0) be a point of order 2 on a nonsingular elliptic curve E 
defined by y^ — f{x) = a;^ _|_ + B over the field ¥q, with quadratic character x- 
Let n be the number of roots of f{x) in ¥q, and for n — i, let xi and X2 denote the 
other two roots. 

For q = 3 mod 4; 

(1) If xif'i^o)) = 1 then E has a point of order 4- 

(2) Otherwise, E has a point of order 4 if and only if n = 3 and xif'i^i)) = 1- 
For q = 1 mod 4; 

(1) If n — 1 then E has a point of order 4 if and only if xif i^o)) — 1- 

(2) Otherwise, if xif i^o)) = 1 {resp., xifi^o)) = —1) then E has a point of 
order 4 if and only if xi^Q ~ ^i) — 1 [resp., xi^i ~ ^2) = !)• 
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Proof. Note that f{xi) — implies f'{xi) ^ 0, since E is nonsingular, hence 
x{f'(xi)) = ±1. Let E denote the quadratic twist of E over F^. By Lemma [1] 
each root Xi of f{x) for which xif'i^i)) = 1 yields 4 points of order 4 (two pairs 
of inverses), either all on E, all on E, or split 2-2 between them. Recall that 
^E = q + 1 — t and = q + 1 + t, where t is the trace of Frobenius, so 4|#i? if 
and only if 4|#i^, and for q = 3 mod 4, 8\#E if and only if 8\#E. 
We first consider q= 3 mod 4. 

Suppose xif'i^o)) = 1- If = 1 then i? and E each have 2 points of order 4. If 
n — 3 then at least one of E and £' has order 8, but if one docs, then so must the 
other, and again E has a point of order 4. 

Suppose xif'i^o)) = If ?i = 1 then E cannot have a point of order 4, so 
assume n = 3. By Lemnia[2l for q = 3 mod 4 we have x(/'(2;i)) = xif'i^^)), and 
if their common value is -1 then E cannot have a point of order 4. If it is 1 then at 
least one of #£' or ^E is divisible by 8, but then they both are and both contain 
a point of order 4. 

We now consider q = 1 mod 4. 

If n = 1 then E can have a point of order 4 if and only if xifi^o)) — Ij as above. 
Now assume n = 3. It follows from Theorem 4.2 of lOJ that E has a point of order 
4 if and only if at least two of xq — xi, xi — X2, and X2 — are squares in Fg. 
We have /'(xq) = (xq — xi){xq — X2), so if xif'i^o)) = 1 then it suffices to check 
x{xa — xi), and if xif'i^o)) — —1 then it suffices to check xi^i ~ ^2)- 

□ 

Lemma 2. Let f{x) be a monic cubic polynomial with distinct roots Xo,Xi,X2 in 
¥q, with q odd. We have 

x(-l)x(/'(a^o))x(/'(x2))x(/'(^2)) = 1. 

In particular, the number of squares in the set {/'(xq), f'{xi), /'(X2)} is even when 
q = 1 mod 4 and odd when q = 3 mod 4. 

Proof. Recall that for a monic / of degree n = 3, the discriminant of / is given by 

D{f) - (-ir("-i)/2i?(/,/') = -R{f,f'), 

where R{f, /') is the resultant. Since / is monic, we have R{f, /') Y[ f'i^i)^ thus 

D{f) = -/'(xo)/'(xi)/'(x2). 

The roots of / are distinct, so D{f) ^ 0. By the Stickelberger-Swan Theorem 
(Corollary 1 in [19 ), D{f) must be a square in Fg, since / is degree 3 and has 3 
irreducible factors. The lemma then follows, since x{D{f)) — 1. □ 

As a final remark, we note that when (1) fails to hold in Proposition!!] it is quite 
likely that E has trivial 4-torsion. On average, this probability is about 90% (this 
can be computed precisely, see [HIS]). As a practical optimization, when seeking a 
point of order 47V, if condition (1) fails we may simply discard the curve and test 
another. When q = 3 mod 4 this reduces to a test for quadratic residuacity in ¥q, 
and we expect two tests of curves generated with Xi{2N) will suffice to produce a 
curve with a point of order 4iV. 
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N F{r, s) 

8 rs-2r + l 

9 r - + s - 1 

10 rs^ — 3rs + r + 

11 — rs^ + 3rs^ — 4rs + s 

12 r^s - Sr^ + r.s + 3r - - 1 

13 - r'^s'^ + 5r^s^ - 9r^s^ + 4r^s - 2r'^ - rs^ + Grs^ - 3rs + r - 

14 r^s^ - 5r^s^ + 6r^s - + rs^ - 3rs^ + 6rs^ - 7rs + r + s 

15 - r^s^ + Tr^s"* - ISr^s^ + IQr^s^ - lOr^s - rs^ + 4rs^ - Srs^ + 5rs - 
+ s'* - + - s 

16 r^s^ - 4r^s + 2r^ + 3r2s^ + 2r2s - 2r^ - rs^ + 4rs^ - lOrs^ + 6rs^ - 3rs 

+ r + s^ 

17 - r'*s<^ + Qr^s^ - 31r'^s'^ + SOr^'s^ _ 39r4s2 + lOr^s - 3r4 - r^s^ + SrS^s 
+ 12r^s^ - 46r^s^ + 54r^s2 _ i^,^3g ^ 3^3 _ ^2^6 _ 3^2 ^5 ^ g^2g4 ^ ^2^3 

- 21r2,s2 + Gr^s - + rs'^ - 3rs*5 + 6rs^ - lOrs^ + llrs^ - 

18 r^s^ _ gj,4^2 _^ gj,4_g _ ^4 _^ j,3_g5 _ 7j,3_g4 20r^ - IQf^s^ - 8r^s + + r^s^ 

- llr^s^ + 28r^s^ + rs^ - 5rs^ - 8rs^ + + + 

19 r*^ - + llrS*5 - 48r5s5 + lOSr^s^ - 121r5s3 _^ gOr^s^ - 20r^s - 

- 2r^s^ + 12r''s'^ - Qr^s^ - QQr'^s^ + 144r'^s3 _ loSr^gS _^ 35^,4^ _ 3^3^? 

+ 3r^s^ + 21r^s^ - 30r^s^ - 41r^s3 + Slr^s^ - 21r^s + r^s^ - Gr^s^ + 21r2s'^ 

- SOr^s*^ + eer^fiS _ Slr^s* + 25r^s^ - 18r2s2 + Tr^s + Zrs^ - Ihrs^ + lOrs* 

- 6rs^ + 3rs^ — rs + 

20 r^s^ - 5r^s + Sr^ + Sr^s^ - lOr* - r^s^ + 9r^s^ - 35r^s^ + 70r^s* - Shr^s^ 
+ blr^s^ - 9r^s + lOr^ + lOr^s^ - 'ihr'^s^ + SOr^s^ - SOr^s^ + lOr^s - 5r^ 

- rs' + 2>rs^ - 6rs^ + lOrs"* - Ibrs^ + ISrs^ - 3rs + r - 

21 r'^ - r^s^ + IZr^s^ - mr^s^ + lQ2r^s^ - ZQQr^s'^ + 2&lr^s^ - llQr^s^ 
+ 21r^s - 4r^ - r^s^ + lOr^s^ - 45r''s'' + 141r4s6 - 3A5r^s^ + 576r^s* 

- 551r^s3 ^ 27Sr'^s^ - 49r'*s + Gr" - r^s^O + lOr^s^ - blr^s^ + IbQr^s"^ 

- 316r^s*^ + 450r^s^ - 551r^s^ + 489r^s^ - 247r^s2 + 42r^s - 4r^ + Sr'^s^ 

- Slr'^s'^ + 109r2s6 _ I72r^s^ + 203r^s^ - ISlr^s^ + 97r^s^ - 14r^s + 

+ 2rs^ - Wrs' + Srs^ + 2rs^ - 13rs'* + 19rs3 - 14rs2 + j-g + gS _ ^7 _|_ ^6 _ ^5 
+ - + 

22 r^fiS _ 9^6 g4 _^ 28r6s3 - 'ihr^ s"^ + ISr^s - + r' 3^ - 12rS^ + 59r5s*5 

- 148r5s^ + 205rS'* - V&^r' + 133rS2 _ 49^,5^ _^ 3^5 ^ ^4^8 _ g^4g7 

- 8r^s'5 + 118r^s5 - 260r*s* + 249r^s3 - 164r*s2 ^ 53^4^ _ 3^4 ^ ^3^8 

- 30r3s'^ + Z\r^s,^ + 70r3s^ - lOGr^s^ + 'mr^s'^ - Z<dr^s + + r^s^ 

+ er^s^ - 7r2sS - 2W'^s^ + Sr^s" + Wr^s^ - ISr^s^ + 7r2s - rs^ + 3rs8 

- 8rs'^ + 21rs'5 - ISrs^ + lOrs^ - 6rs^ + Zrs^ -rs-s' 

23 - r^s^ + ISr^s^ - 94r8s'' + Sigr^s^ - %mr^s^ + 756r^s4 - 520r^s3 

+ 189r8s2 _ 35^8^ _ _ 4^7g9 _^ 39^.7^,8 _ i20r7s'^ + 28r^s6 + 597r^s5 

- 1341rV + 1256r^s3 - 525r''s2 + lOSr^.s + - lOrS^ + AS^r^ + 24rSs^ 

- 357r''s6 + 324r6s5 + hlQr^s^ - 1130r6s3 + 576rS2 - 126r6s + r^s^^ - I4r5si 
+ 93r5s" - 370r5sio + 970rS9 - 1827r5s8 + 2553r\s^ - 22'd%r's^ + lOgSr^s^ 

- 480r5s'' + eSGr^s^ - 369rS2 + 84r'5s + r^s^^ _ 21r4s" + leSr^s^o - GSOr^s^ 
+ 1530r4s« - 2562r4s^ + 2957r4s6 - 2046r4s-^ + TSOr^s* - Mhr^s^ + 171r^s2 

- 36r'*s + r\si2 - ISr^s" + %%r^s^^ - Mr^:^^ - 45r3s^ + 402r3s'' 

- 833r3s*^' + 837r-\s''^ - SSlr^s'' + \A5,r^s^ - A&r^s'^ + 9t^s + r^s^^ - gr^s" 
+ 13r2si° - r^s^ - 2Ar^s^ + 2Sr^s'^ + A2r^s^ - UQr^s"^ + SGr^s^ - 21r2s-^ 
+ er^s^ - r^s + rs^^ _ 3^,^11 _^ 5^,^10 _ ;^q^^9 ^ ^g^^s _ 21^57 _^ 21rs^ - 



Table 4. Raw form of Xi(Af) : F{r, s) = 0. 
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J/ \ 

d{Co) 


d{Ci) 


t{Co) 


t{Ci) 


^max 


^(Co,Ci) 


10 





4 





1 


1 


2 


10 


11 


1 


2 


2 


5 


4 


2 


4 


12 





2 





6 


1 


2 


13 


13 


2 


3 


2 


11 


6 


2 


13 


14 


1 


2 


2 


10 


4 


2 


11 


15 


1 


3 


2 


15 


5 


3 


18 


16 


2 


3 


2 


13 


6 


5 


23 


17 


5 


5 


4 


28 


12 


5 


23 


18 


2 


4 


2 


19 


6 


5 


24 


19 


7 


6 


5 


39 


18 


4 


23 


20 


3 


5 


3 


28 


6 


4 


23 


21 


5 


6 


4 


55 


11 


4 


18 


22 


6 


6 


4 


50 


17 


7 


40 


23 


12 


9 


7 


87 


38 


7 


25 


24 


5 


6 


5 


41 


20 


6 


25 


25 


12 


10 


8 


114 


46 


6 


20 


26 


10 


8 


7 


82 


27 


5 


32 


27 


13 


11 


8 


135 


52 


4 


19 


28 


10 


10 


7 


115 


30 


2 


16 


29 


22 


14 


11 


214 


88 


8 


32 


30 


9 


10 


8 


109 


46 


7 


23 


31 


26 


16 


13 


279 


124 


6 


23 


32 


17 


13 


10 


190 


78 


7 


19 


33 


21 


16 


12 


319 


109 


6 


29 


34 


21 


14 


11 


235 


88 


7 


22 


35 


25 


19 


15 


438 


142 


4 


19 


36 


17 


14 


11 


224 


94 


7 


23 


37 


40 


23 


18 


582 


225 


4 


19 


38 


28 


18 


14 


383 


140 


6 


27 


39 


33 


22 


17 


586 


212 


4 


20 


40 


25 


19 


15 


412 


171 


5 


22 


A 1 

41 


oi 


zo 


oo 


( U 




Q 
O 




42 


25 


20 


15 


442 


165 


8 


27 


43 


57 


31 


24 


1065 


408 


6 


23 


44 


36 


24 


19 


654 


208 


3 


21 


45 


41 


29 


23 


960 


368 


4 


19 


46 


45 


26 


21 


791 


285 


6 


23 


47 


70 


37 


29 


1526 


1768 


6 


33 


48 


37 


26 


19 


773 


257 


7 


23 


49 


69 


39 


31 


1791 


900 


6 


37 


50 


48 


30 


23 


1040 


391 


8 


42 



Table 5. Search algorithm statistics for Xi{N). 



Co and Ci are (respectively) the raw and optimized forms of Xi{N). The column d(Ci) 
denotes the minimum of the degree of d in x or y, and t{Ci) denotes the number of terms. 
The column l{Co, Ci) gives the length of the path traveled by the algorithm of Section |3] 
to reach Ci from Co (typically not a shortest path), and fcmax is the maximum value of k 
prior to reaching Ci. 
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N f{x,y) 

11 y'^ + {x^ + l)y + X 

13 y"^ + {x^ + x"^ + l)y ~ x^ ~ X 

14 + {x^ + x)y + x 

15 y^ + {x^ + X + l)y + x^ 

16 y'^ + {x^ + x^ — X + l)y + x"^ 

17 y* + {x^ + x'^ ~ X + 2)y^ + {x^ - 3x + l)y^ - {x* + 2x)y + x^ + x'^ 

18 y^ + ix^ -2x^ + 3x+l)y + 2x 

19 y^ - {x^ + 2)y-* - {2x^ + 2x^ + 2x - l)y^ + {x^ + Sx"* + 7x^ + 6x^ + 2x)y'^ 

-{x^ + 2x^ + 4a;^ + ?,x'^)y + x^ + x^ 

20 y3 _^ {^^2 _^ 3^y2 _^ (^3 _^ _^ 2 

21 y4 + (3a;2 + 1)^^ + [x^ + x'' + 2j;2 + 2x)y'^ + {2x^ + x^ + x)y + x^ 

22 y-* + [x^ + 2x'^ + X + 2)y^ + {x^ + a;^ + 2x^ + 2x'^ + l)y'^ 

+ {x^ — x'^ — 2x^ — x"^ — x)y — x"^ — x^ 

23 y^ + {x^ - a;^ + + Ax'^ + 3)y^ + (x^ + Sx^ + + Sx^ + Tx^ - 4a: + 3)y'5 

+ {2x'^ + ix^ -x^ - 2x^ - - 8x + l)y'' 
+ - 4a;6 - hx^ - - ^x^ - 2x'^ - 3x)y^ 

-{3x'^ - 5a;4 - - Sx^ - 2x)y'^ + (3a;^ + Ax'^ + x)y ~ x'^{x + if 

24 yS + + Ax^ + 3x'^ - x - 2)y^ - {2x* + 8x^ + Ix^ - l)y^ 

-{2x^ + 4a;4 - Sx^ - hx'^ - x)y^ + {2x^ + Sa;* + 2a;^)y + a;^ + a;^ 

25 yS + (4a;2 + 7a; - 4)y^ - (a;^ - .x* - Ua;^ - 4a;2 + 24a; - <a)y^ 

-(a;'^ + Ax^ ~ 3a;^ - ISa;^ + ISx^ + 33a;2 - 30a; + A)y^ 

-{x^ + 2a;'^ - 8a;^ - 14a;-'5 + 24x'' + 17a;3 - 41a;2 + 16a; - l)y4 

+ (.a;^ + 6.a;'^ + 3a;*^ - 20a;-'^ - 3.x'' + 28x3 - lOx^ + 3x)y3 

-(3x'^ + Ox*^ - 3x-'^ - 13x'' + llx^ - 3x2)y2 + (3x6 ^ 4^5 _ 4^4 _^ ^3)^ _ ^5 

26 yS + (3x2 + 42; _ 2)y5 + (3x4 ^ ^Qx^ - 9x + l)y4 

+ (x6 + 7x5 _^ g^4 _ ;^4^3 _ Y\x'^ + 6x)y3 
+ (x^ + 4x6 - x^ - 13x4 + 2x3 ^ ;lox2 - x)y2 
-(x^ - 7x4 _ 4^3 _,_ 2x'^)y - x4 - x^ 

27 yS + (3x2 + 6x - 3)y^ - (3x5 _ 2^83,3 _ 93.2 _^ ^^g^; - 3)y6 

-(x^ + 8x'^ + 13x6 - 21x5 - 48x4 + 20x3 + 42x2 - 18x + l)y5 

-(xi° + 6x9 ^ ^2x® - 14x^ - 72x6 _ 27x5 + 93x4 ^ 33^3 _ 45^2 ^ g^)y4 

+(xi° + llx^ + 40x** + 36x'^ - 69x6 _ ^^^x^ + 33x4 ^ 54^3 _ 152:2)2^3 

-(4x9 + 30x** + 63x^ + 10x6 - 69x5 _ 24a;4 ^ V^x^)y'^ 

+ (6x^ + 27x'^ + 27x6 - 6x5 „ i2a;4)y _ Zx'^ _ 6x6 _ 33,5 

28 y'^ + 3xy6 + (x^ + 3x4 ^ 5^3 _^ 9^.2 _^ 2x)y5 - (2x5 - 6x3 ^ 2x2 + 2x)y4 

+ (3x6 _^ ;l6x5 + 18x4 - 2x2)y3 + (x^ - 2x6 „ 20x5 - 28x4 - 12x3 - 2x2)y2 
-(2x'^ + 3x6 _ 52,5 _ iQ^^ _ 5^3 _ ^2^y + _^ 2x6 ^5 

29 y" + (2x3 _^ 5^,2 + 5^; _ 3)^10 (^6 ^_ ^g^i _^ ^^^s _ _ i2x + ■ ■ ■ 

30 y* - (2x3 + 4x2 + X + 5)^"^ + (x6 + 4x5 ^_ g^4 _,_ 9^3 _,_ ^^4^2 ^_ ^o)y6 

-(x^ + 4x6 + 9x5 _,_ iQ^i ^ 4^,3 i^^2 _ iQ^y5 

+ (x^ + 4x^ + 4x6 - 5x4 _ 20a;3 + 5a;2 _ 20x + 5)y4 
+ (3x'^ + 11x6 + 15x5 + 9x4 _,_ -^3^3 _ 9^,2 ^4^ _ i^yS 

+ (3x6 _^ 9^5 _^ ^42-4 ^ 2x3 _^ i^^2 _ 32.)y2 ^ (^5 + a;4 _^ 42;3 _ 32-2 )y _ 2.3 



Table 6. Optimized form of Xi{N) : f{x,y) = 0. 

The polynomial for TV = 29 is not displayed in full. Full polynomials for A'^ < 50 are available at 

|http : //math.mit ■ edu/^drew. 
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N 

6 r = X, s = 1 

7 r = X, s — X 

8 r = l/{2-x), s = x 

9 r = x'^ — X + 1, s = x 

10 r = -a;2/(x2 - 3a; + 1), s = x 

11 r = 1 + xy, s = \ — X 

12 r = (2x2 - 2a; + 1) /x, s = (Sx^ -3x+l)/x'^ 

13 r = 1 — xy, s = 1 — xy/{y + 1) 

14 r= l-(x + y)/((y+l)(x + 2/ + l)), s = (1 - x)/(2/ + 1) 

15 r = 1 + (xy + y'^)/{x^ + x'^y + x'^), s = 1 + y/{x'^ + x) 

16 r = (x^ — xy + + jy)/(x^ + X — y — 1), s = (x — y)/(x + 1) 

17 r = (x^ +x - j/)/(x2 +xy + x - - y), s = (x + l)/(x + y + 1) 

18 r= (x^ -xy-3x + l)/((x- l)2(xy + l)), 
s = x^ — 2x — y)/(x2 — xy — 3x — y^ — 2y) 

19 r = 1 + x(x + y)(y - l)/((x + l)(x2 - xy + 2x - y^ + y)), 
s = 1 + x(y - l)/((x + l)(x - y + 1)) 

20 r = 1 + (x^ + xy + x)/((x - l)'^{x'^ -x + y + 1)), 
s = 1 + (x^ + y + l)/((x - l)(x2 - X + y + 2)) 

21 r = 1 + (y2 + y)(xy + y + l)/((xy + l)(xy - y^ + 1)), 
s = l + (y2 + y)/(xy + l) 

22 r = (x^y + x2 +xy + y)/(x3 + 2x2 + y), s = {xy + y)/{x'^ + y) 

23 r = (x2+x + y + l)/(x2-xy), s=(x + y + l)/x 

24 r = (x^ + X - y + l)/(x2 + xy - y2 + y), s = (x + l)/(x + y) 

25 r = (x^ +xy + y2 - y)/(x2 +x + y - 1), s = (x + y)/(x + 1) 

26 r = (x^y + 3x2y - x^ + xy^)/{{x + l){x'^y + x^ + 3xy + y^)), 
s = (xy - x)/(xy + y) 

27 r = (— x^ — x^ — X — y)/(x2y + xy — X — y), s = {—x^ — x — y)/{xy — x — y) 

28 r = 1 + (xy + y)/((y - l)(xy - x + 2y - 1)), 
s = 1 - (xy + y)/((y - l)(x - y + 1)) 

29 r = (— x"^ — x^ — X — y)/{x^y + xy — x — y) 
s = 1 - (x^ + xy)/{xy - x-y) 

30 r = (x^y + X + y)/ (x^y — xy + x), 
s = (x^y + xy + X + y)/ (x^y + x) 



Table 7. Birational maps for Xi{N) from f{x,y) = to F{r,s) = 0. 



